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Developmental Partial Differential Equations 


Nastassia Pouradier Duteil 1 , Francesco Rossi 2 , Ugo Boscain 3 , Benedetto Piccoli 1 


Abstract — In this paper, we introduce the concept of Develop¬ 
mental Partial Differential Equation (DPDE), which consists of a 
Partial Differential Equation (PDE) on a time-varying manifold 
with complete coupling between the PDE and the manifold’s 
evolution. In other words, the manifold’s evolution depends on 
the solution to the PDE, and vice versa the differential operator 
of the PDE depends on the manifold’s geometry. DPDE is used 
to study a diffusion equation with source on a growing surface 
whose growth depends on the intensity of the diffused quantity. 
The surface may, for instance, represent the membrane of an 
egg chamber and the diffused quantity a protein activating a 
signaling pathway leading to growth. Our main objective is to 
show controllability of the surface shape using a fixed source 
with variable intensity for the diffusion. More specifically, we 
look for a control driving a symmetric manifold shape to any 
other symmetric shape in a given time interval. For the diffusion 
we take directly the Laplace-Beltrami operator of the surface, 
while the surface growth is assumed to be equal to the value of 
the diffused quantity. We introduce a theoretical framework, 
provide approximate controllability and show numerical results. 
Future applications include a specific model for the oogenesis 
of Drosophila melanogaster. 

INTRODUCTION 

Many mathematical models aim to reproduce biological 
mechanisms, including those involving growth of living 
organisms (see [3], [7], [16], [18]). This concerns the field 
of Developmental Biology, which is devoted to the study of 
growth and development of organisms, as well as the genetic 
control of cell growth, differentiation and morphogenesis 
(see [21], [25]). 

We focus on providing a mathematical framework for 
understanding growth of organisms induced by signaling 
pathways. Many approaches for modeling biological growth 
were proposed in the literature, see for instance [10] and 
references therein. In that case, a tissue is either regarded 
as a continuum or as a collection of cells. The latter 
microscopic approach is based on discrete models such as 
cellular automata. On the other hand, Partial Differential 
Equations are natural for the former macroscopic point of 
view, where the dependent variable usually represents mass 
concentration. In most models the growth of mass is then 
assigned as an internal mechanism, e.g. proportional to the 
mass concentration, or as an external one. The point of view 
is similar to that of elasticity, but the equations are modified 
to include growth (so violating conservation of mass). 
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There exist abundant investigations on Turing Patterns 
(TP), which are generated by two chemicals due to instabili¬ 
ties caused by diffusion (see [23], [24]). Papers also investi¬ 
gated effects on stability, geometry and growth in TP due to 
growth rate, curvature and other characteristics of chemical 
interactions ([2], [4], [6], [8], [14], [15], [17], [20]). However, 
to our knowledge, there is no clear experimental evidence of 
model organisms showing patterns via mechanisms predicted 
by the corresponding mathematical model. 

In our approach, we aim at having a faithful representation 
of the fact that growth is mostly regulated by signaling 
pathways. In order to achieve this, we design a mathematical 
framework based on two main ingredients. The first is a 
smooth topological manifold M (representing for instance 
an egg chamber, a tissue, or even a single cell) which 
evolves in time. The second is a quantity s evolving on the 
manifold according to a specific PDE. More precisely, we 
assume that the manifold’s change (usually growth) depends 
directly on the quantity modeled by the PDE and, in turn, the 
PDE operator depends on the (changing in time) manifold’s 
geometry. In order to define the model in a mathematically 
sound way, we consider a manifold M embedded in a 
ambient Euclidean space M n and parameterized (possibly 
locally) by variables y E M m via a map ^(t^y). Moreover, 
the evolution of M is given by a vector field u(s, x), x E M n , 
depending on 5. To mimic the genetic mechanisms regulating 
growth, we introduce a control term which sets the intensity 
of a source for the quantity s. Finally the coupled system 
reads: 

X M s = S(u), (1) 

where Xm is a differential operator defined on the manifold 
M, S a source term and u the control. 

PDEs on manifolds form a vast subject, with a wide and 
varied literature. Previous works focusing on the control of 
PDEs on manifolds include [5], [11]. Other works have dealt 
with many different PDEs on manifolds. For instance, [22] 
focused on harmonic maps on Riemannian surfaces, while 
the Cauchy problem for wave maps is treated in [19]. A 
general approach for hypoelliptic Laplacian on unimodular 
Lie groups can be found in [1]. Regarding controlled evo¬ 
lution on manifolds, the Klein-Gordon equation on a 3-D 
compact manifold is considered in [13]. All these works 
focus on PDEs on manifolds (and possibly their control) but, 
to our knowledge, our contribution is the first that deals with 
a completely coupled system manifold-PDE of the type ([]]). 
Because of its biological meaning, we call such system a 
Developmental Partial Differential Equation (DPDE). 



The specific application we have in mind is that of under¬ 
standing how the Gurken protein contributes to the formation 
of various Drosophila eggshell structures. During oogenesis, 
Gurken is secreted near the oocyte nucleus, and then diffuses 
and is integrated by the epidermal growth factor receptor 
(EGFR), triggering a signaling pathway. Therefore, modeling 
evolution of the Gurken concentration during Drosophila 
oogenesis is more complex than solving a PDE on a manifold 
due to the fact that the egg chamber grows throughout 
the process. Leveraging on natural symmetries of the egg 
chamber we choose as M a one-dimensional symmetric 
manifold embedded in M 2 and initially equal to S 1 . For 
the operator Xm, we simply choose the Laplace-Beltrami 
operator on M . Finally, the resulting DPDE is given by ([9]). 
Our main aim is to show controllability in terms of the 
possible shapes reachable from S 1 regulating one or more 
sources. We show how to adapt the approach of Laroche, 
Martin and Rouchon [12], proving flatness of the heat 
equation, to our setting and then provide numerical studies. 

The paper is organized as follows. We begin by introduc¬ 
ing the biological context that motivates our work, that is 
diffusion of Gurken in the growing Drosophila egg cham¬ 
ber during oogenesis. We then establish the mathematical 
framework, starting with a general setting involving multiple 
reaction-diffusion equations on a dynamically evolving man¬ 
ifold. We simplify the problem by considering the diffusion 
of just one quantity (the growth signal) and its direct impact 
on the manifold’s evolution. Finally, we show numerical 
simulations depicting various possible shapes that can be 
obtained by controlling one or two sources for the signal. 

I. BIOLOGICAL CONTEXT 

Drosophila melanogaster is a commonly used model or¬ 
ganism to study cell signaling, tissue patterning and mor¬ 
phogenesis. Tissue patterning and cell fate determination are 
guided by a handful of cell signaling pathways. For instance, 
the epidermal growth factor receptor (EGFR) signaling 
pathway controls many cell processes, including apoptosis 
(cell death) and cell migration. In particular for Drosophila 
oogenesis, the TGF-alpha like ligand Gurken is secreted from 
near the oocyte nucleus, diffuses in the perivitelline space 
surrounding the oocyte, and signals through EGFR in the 
overlaying follicle cells. This will eventually gives rise to 
cell differentiation, forming various eggshell structures. 

These processes occur over the course of about 27 hours, 
during which the Drosophila egg chamber undergoes mor¬ 
phological changes, including the follicle cells’ gradual 
movement over the oocyte (see Figure 1.). Thus, differ¬ 
ent cells are dynamically exposed to Gurken over time. 
Mathematically, this is equivalent to considering the oocyte 
nucleus as a moving source of Gurken. Moreover, during 
oogenesis the egg chamber grows and becomes wider and 
more elongated. 

The diffusion of Gurken has been previously modeled to 
determine its impact on local genes (see [9], [26]). These 
existing models consider diffusion at steady-state from a 
fixed nucleus, thus not accounting for the evolution of the 


egg chamber’s shape. The aim of this work is to propose a 
mathematical setting that will allow to model the growth of 
the egg chamber as oogenesis progresses. 



(a) Schematic of the eggshell growth during oogenesis from Stage 6 to Stage 
11. The resulting egg with protruding dorsal appendages is shown on the right. 



(b) The Gurken ligand diffuses in the oocyte and binds to EGFR, triggering a 
signaling cascade. The level of EGFR signaling can be assayed with dpERK 
antibody staining. The non-homogeneous activation of EGFR is responsible 
for the formation of various morphological features on Drosophila species’ 
eggshells, such as the dorsal appendages or the dorsal ridge. 


Fig. 1: From [26] J.J. Zartman, L.S. Cheung, M. Niepielko, 
C. Bonini, B. Haley, N. Yakoby and S.Y. Shvartsman, 
Pattern formation by a moving morphogen source, Phys 
Biol 8 (2011) 045003. ©IOP Publishing. Reproduced with 
permission. All rights reserved. 


II. MATHEMATICAL FRAMEWORK 
A. GENERAL SETTING 

One way to model the growth of a manifold (for instance 
representing the Drosophila egg chamber) is to consider that 
it is provoked by a signal diffused from a source. We assume 
that the manifold is embedded in M 3 and parametrized 
by the two-dimensional sphere S 2 . We indicate by r the 
manifold radius, which is a function of time t (and of the 
S 2 coordinates) and evolves depending on the signal s : 

d t r = F(s). (2) 

The signal s also evolves with time following the reaction- 
diffusion equation: 


d t s = A r s + R s {s), (3) 

where A r denotes the Laplace-Beltrami operator on the 
surface given by r. Other signals may evolve on the manifold 
according to reaction-diffusion equations but not inducing 
growth. They can be modeled by: 

dtrrii = A r nrii + Ri{rrii). 


(4) 




B. TOY MODEL 

We simplify the problem described in ([2]), ^ and ([4]) by 
focusing on a one-dimensional manifold and by neglecting 
the proteins rrii that do not contribute to the growth of the cell 
membrane. Given an angle variable 9 E S 1 , we characterize 
the position of the membrane by a function r = r(t,9) 
representing the radius. Moreover, we consider that F = Id, 
so that the signal s = s(t, 9) directly pushes the manifold to 
grow in its radial direction. The dynamics of s is given by 
the heat equation on the manifold. We neglect the reaction 
term, compensating for it by allowing the signal 5 to become 
negative. We introduce a control u , which sets the value of 
s at the point 9 = it. Biologically, this corresponds to the 
point at which the nucleus sends the growing signal to the 
boundary. Hence, the dynamics satisfies: 

d t r = s, 

^ dts — A r s, (5) 

^s(£, 9 = 7r) = u(t). 

We now assume that the initial configurations of both r 
and s are symmetric with respect to 6, i.e. r(0, — 9) = r(0, 9) 
and, similarly, s(0, — 9) = s( 0,0). The simplest example is 
r(0,$) = 1 and s(0,9) = 0, i.e. a round cell and a zero 
signal on it. One can easily prove by that, for any choice of 
the control u(t ), both r and s stay symmetric. Indeed, using 
the explicit expression of the Laplace-Beltrami operator ([9]), 
we prove that (s(t, #), r(t, 9)) and ((s(£, — 0), r(t, —6)) solve 
the same differential system. Since the two couples have 
identical initial conditions, by uniqueness of solution we 
deduce that they are equal and thus symmetric. 

Since s is the solution of a heat equation, it is a C°° 
function far from 9 = tt for all time. As a consequence, 
symmetry also implies dos(t,Q) = 0 for all t. Hence, we 
reduce our study to the half-circle 6 E [0, tt] and consider 
the following dynamics: 

' d t r = s, 

< dtS = ArS ' ( 6 ) 

s(t, 9 = 7r) = u(t), 

< d#s(t , 9 = 0) = 0. 

We now study the Riemannian structure on the cell in¬ 
duced by a shape r. As already stated, s is C 00 except in 
0, since its value there depends on u(t). Assuming that the 
choice of u implies that s is C 00 at 0 too, we have that r is 
a C°° function too. Consider the coordinate 9 on the circle, 
and observe that a displacement do on the coordinate induces 
a displacement in the r variable that can be estimated by 
y/r 2 + rfid#, where tq = dor is the derivative of r with 
respect to 9. The estimate is due to a simple geometric 
first-order estimate of the length of the curve r(9). As a 
consequence, one can define the following metric on S' 1 : 

go is bilinear and satisfies: go(do-,do) = r 2 (9) -\-rg(9). (7) 

This uniquely defines the metric on S 1 . It is also clear that 
the inverse of the metric satisfies g e (d9,d9) = r 2(o)+ r 2 (d) • 


Such an operator is never zero since the radius is supposed 
to be positive for all 9. Then, a direct computation gives the 
explicit expression of the Laplace-Beltrami operator A r . We 
have: 


A r s = 




y'r 2 +^ d> (yV + rJ 


r z + r; 


:d 2 e s- 


rro 


r -^a ts . 


(r 2 + rg) 2 


( 8 ) 


Hence the system we want to study is the following: 


d t r = s, 
d t s 




rro+rodgr 

( r 2 _|_ r 2)2 


dos , 


s(£, 9 = tt) = u(t ), 
d e s(t,9 = 0) = 0. 


(9) 


We want to prove controllability for system in a 
specific case, that is to find a control u that drives a 
(symmetric) cell shape to another (symmetric) cell shape in 
a given time interval [0,T], together with having a signal 
s that is zero at the initial and final times. In mathematical 
terms, we consider initial and final configurations ro, r\ and 
a time T > 0. We want to find a control u : [0; T] M 
such that the unique solution of with r(t = 0) = ro and 
s(t = 0) = 0 satisfies r(t = T) = n and s(t = T) = 0. 
This goal is called exact controllability. It is known that this 
goal is impossible to be achieved in general, since we already 
know that some configurations (for instance non-smooth final 
configurations) cannot be reached with a heat equation. 

Hence, we instead aim to prove approximate control¬ 
lability , defined as follows: considering initial and final 
configurations ro, r\ and a time T > 0, for every e > 0, 
we want to find a control u : [0;T] —>> M such that the 
unique solution of ([9]) with r(t = 0) = rO and s(t = 0) = 0 
satisfies || r(t = T) — ri \\ L 2 < e and || s(t = T) \\ L 2 < e. 

It was shown in [12] that the 1-D generalized heat equation 


' dt4> = f(6)dg<t> + g{0)de4> + 

< cj){t, 9 = 7r) = u(t), (10) 

do<p(t , 9 = 0) = 0 


is approximately controllable where / > 0, g and h are an¬ 
alytic functions. Moreover, [12] proves a stronger condition: 
(approximate) motion planning or (approximate) tracking, 
defined as follows. Given a reference trajectory, we want 
to find a control such that the solution of the system © 
stays close to the reference trajectory for each time. In 
mathematical terms, one has the following result. 

Theorem 2.1: Consider a time horizon [0, T] and a smooth 
trajectory / : [0, T] —)> L 2 (0,7r). For every e > 0, there exists 
u : [0, T] such that the solution of (TT)] ) with initial data 
/(0) satisfies || f(t) — f(t )\\ L 2 < e for all time t E [0, T\. 

We use this result to prove approximate controllability of 
Moreover, we will show a stronger condition, that is 

















approximate tracking of the r variable, together with the 
condition || s(t = 0)|| L 2 < e and \\s(t = T )\\ L 2 < e. 
Since we need analytic coefficients for the second equation 
of ([9]), we need a reference trajectory that is analytic for all 
t, i.e. r : [0; T] —)• C w (0,7r), together with smoothness with 
respect to t. We can prove the following main theorem. 

Theorem 2.2: Let f : [0, T] C w { 0,7r) be a reference 

trajectory. Then for all e > 0, there exists a control u : 
[0, T] M such that the unique solution of © with r(t = 
0) = r(t = 0) and s(£ = 0) = 0 satisfies ||r(£) — f (£)|| L 2 < e 
for all t G [0, T]. 


III. EQUILIBRIA 


We look for equilibria of the form: u(t ) = u e and s(£, 0) = 
s e ($), that solves the system: 


d t r e = s e , 

f)2 „ _ r e d e r e +d e r e d%r e 

U 6 b e — (r e ) 2 +(d e r e ) 2 ° ebe 

S e (0 = 7r) = U e , 

%s e (0 = 0) = 0. 


( 11 ) 


From ([9]), we deduce that for all 6 , r e (£, #) is a linear function 
of s e ($): 

r e M) = s e (0)t + r o (0), (12) 


where ro(0) := r e (0,6 > ). One obvious possible equilibrium 
is obtained when there is no control, i.e. for a zero signal 
(since s e then solves a Laplace equation with the boundary 
condition s e (0 = 7r) = 0). One gets: 

u e = 0 , 

< s e (6) = 0 for all 0 G [0,7r], 

j e (t,6) = ro(0) for all t G [0,T], for all 0 G [0,7r]. 

Hence, if s e and u e are at an equilibrium such that u e = 0, 
there is no signal and the radius is constant in time. 

On the other hand, if u e > 0, then s e solves a Laplace- 
type equation with a non-zero Dirichlet boundary condition 
at 6 = 7T, so s e (6) > 0 for all 0 G [0,7r]. Hence r e (t, 0) grows 
linearly with time and does not reach an equilibrium. We 
instead look for an equilibrium in the shape of the membrane, 
by defining p e (t,6) = r ( n °ti ce that this is possible 

since r e ^ 0). Then p e is constant in time if d t p e = 0, which 
gives: 


d t r e (t, 9)r e (t , 7r) - d t r e (t, n)r e (t, 0) = 0. 


Since d t r e (t,6) = s e (6 ), we get: 


dtPe(t,9) = 0 


r e (t,9) 
r e (t, tt) 


s e (0) 

SeM' 


This means that at each time t, the membrane r e is a dilation 
of the signal s e . In particular, at t = 0, ro($) = ^^y«s e (6 > ) 

for all 0. Hence from (T~2| ) we get: r e (t, 6 1 ) = g e (^)(f+ )• 

Since s e (6 ) and r e (t, 6 1 ) are proportional, the second equation 
of ( [TT] ) becomes: 


de s e 


SedpSe + doSedpSe 

(S e ) 2 + (SflSe) 2 


which, after simplification, gives: 

S e dgS e = ( d0S e ) 2 . 

One solution to this nonlinear differential equation is the 
constant signal s e (6) = u e , where 5 e satisfies both the 
Neumann and Dirichlet boundary conditions prescribed in 



We relax our conditions and look for a solution s e that 
satisfies s e (7r) = u e but not dos e ( 0) = 0. In particular, if we 
suppose that ds e (0 ) ^ 0 for all 0 G [0,T], we can write: 

dffS e _ d 0 s e 
d 0 s e s e 

Then ^(ln(^s e )) = 9 0 (ln(s e )), so we get: s e {6) = 
u e e x ^ e ~ n \ where A is a constant. Notice that then we 
can bring d 0 s e ( 0) = u e \e~ X7r arbitrarily close to zero by 
choosing A, so we partially recover the original Neumann 
boundary condition. 


IV. SIMULATIONS 

We simulate diffusion of the signal by discretizing the 
second equation of system using Finite Differences, 
supplemented by a Neumann boundary condition at angle 
6 = 0 (d 0 s(t, 0) = 0) and a Dirichlet boundary condition 
at angle 9 = tt (s(£,7t) = u(t)). Then the radius of the 
manifold at each time-step is obtained by simple integration 
of the signal. 


A. Comparison of diffusion on a static vs growing manifold 

We run simulations for a constant control u\ = 1, an 
initial signal so(0) = 0 and an initial radius ro(6) = 1 for 
all 0 G [0,7r]. We notice that s reaches an equilibrium after 
time t = 2. After that point, the radius grows in a linear way, 
i.e. p(t ) = const. See Figure [2] 

We then turn our attention to the comparison with the 
case in which we neglect the growth of the manifold (this 
would correspond to an egg chamber of constant size). In 
this case, taking as initial condition a circle, the radius r 
is constant both w.r.t. time and the 6 variable, thus r = 1 
and r 0 = 0. Plugging this information into equation ([9]), the 
Laplace-Beltrami operator reduces to standard diffusion and 
we get the following system: 


' d t r = s, 
d t s = <9|s 
s(t, 6 = 7r) = u(t ), 

^d e s{t,e = 0) =0. 


(13) 


The simulations for a constant control u = 1 are very 
different from those obtained by using the system &■ 
Figure [3] shows the evolution of the signal and the radius 
with constant control for system. The signal s reaches an 
equilibrium s(t,6) = 1, which means that the growth of 
the radius tends to be uniform with respect to the angle 6. 
Therefore, as expected, neglecting the growth of the manifold 
generates uniform growth and, in the biological system, 
would give rise to spherical egg chambers opposed to the 
spheroidal ones observed in nature. 















Fig. 2: Signal s (left) and radius r (right) for a constant 
control u = 1 at times t = 0.1, t = 2 and t = 8. The source 
correponds to the angle 0 = 7r, so in the signal picture it is 
located on the left end of the equator line corresponding to 
coordinates (—1,0). 





Fig. 3: Signal s (left) and radius r (right) for a constant 
control u = 1 at times t = 1, t = 2, t = 5 and t = 8. 


where uj = ^ so that we obtain a complete sinusoidal 
oscillation up to time 5 then the signal is vanishing (which 
coincides with control u 2 depicted in Figure [4]). The final 



Fig. 4: Control functions u\, u 2 and us 


result is a apple shape manifold with pitch located at the 
signal source point, see Figure [5] To better visualize the 
relationship between the signal and the shape we visualized 
the signal on the manifold itself, so for positive values the 
signal will be outside the manifold and inside for negative 
ones. 


t=1.00 t=3.00 




B. Single source 

A source placed on the first axis (at angle 0 = n) allows us 
to control the diameter of the manifold along the same axis. 
In Figure [2] the manifold is stretched along the first axis 
direction at final time, with an emphasis on the left side, 
i.e. r(T, 7r) > r(T, 0). Using the source to impose negative 
values of the signal (which has a mathematical meaning but 
not a biological one), we can control the final shape of the 
manifold to achieve r(T, 7r) < r(T, 0). In order to do that 
we set the control as: 

f 0.5 • sin(ut) t G [0, 5] 

\ 0 te] 5,10] ’ 


Fig. 5: Radius r (in blue) and signal s (plotted as r + s in 
red) for a control u = at times t = 1, 3, 5, and 10. 

Using a single source it is also possible to induce an 
homogeneous growth along all directions, but with time- 
dependent signals. We first give an impulse and then turn 
off the signal. Define the control by: 

„U\ = f °- 2 ■ sin{ut) t e [0,2.5] 

{ > \ 0 te]5,10] ’ U ; 

where ui = ^ L , so that the half sinusoidal oscillation gives 
an always positive signal (this correspond also to the control 
us depicted in Figure]?]). The final shape is close to that of 


u(t) = 





































































































a circle, but with a larger radius than that at initial time (see 
Figure [ 6 ]). 

t=1.00 t=3.00 




t=5.00 t=10.00 








Fig. 7: Radius r (in blue) and signals sl (plotted as r + sl 
in red) and sr (plotted as r + sr in green) for controls 
ul = ur = U 3 at times t = l, t = 3, t = 5 and t = 10. 


Fig. 6 : Radius r (in blue) and signal s (plotted as r + s in 
red) for a control = 1 x 3 at times t = 1, 3, 5 and 10. 


C. Double source 


As observed above, a single static source allows us to 
control the radii r(T, 0) and r(T, 7 r), i.e. the horizontal 
growth. In order to achieve a larger growth along the vertical 
axis, we consider a system with double source: one locate at 
angle 6 = 0 and the second (as before) at angle 0 = tt. We 
obtain the system: 


d t r = S L + S R , 

1 dtsL-tffiffideSL 

rre+r °M r -des R 


r 2 +r; 


7 e *- 

! - 


d t SL = 

dt$R r 2_^_ r 2 ^ r 2 _|_ r 2^2 

sl(*, 6 = n) = u L (t ), 0 = 0) = ixj*(t), 

dos L (t, 6 = 0 ) = 0 , d e s R (t, 6 = tt) = 0 . 


(15) 


If we use the control given by formula (I4| ) for both sources, 
we obtain a final manifold stretched more in the vertical 
direction, i.e. r(T, 7 r/ 2 ) > r(T, 0) = r(T, 7 r), see Figure [ 7 ] 


CONCLUSION 


In this paper we introduced a new mathematical frame¬ 
work, called Developmental Partial Differential Equation 
(DPDE), to model the growth of organisms induced by 
signaling pathways. A DPDE consists of a couple: a time- 
varying manifold and a signal evolving on the manifold. In¬ 
spired by the specific application to Drosophila egg chamber 
development, we consider a completely coupled evolution 
where the manifold’s growth is regulated by the signal 
and the signal diffusion by an operator (Laplace-Beltrami) 
depending on the manifold geometry. 

We provide controllability results using flatness of the heat 
equation and show simulations of resulting manifold shapes. 


Future work will include: moving sources, general differ¬ 
ential operators, higher dimensions manifolds and explicit 
expressions of controls for motion planning. 
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